Inadequacy of fluvial energetics for describing gravity current autosuspension

Gravity currents, such as sediment-laden turbidity currents, are ubiquitous natural flows that are driven by a density difference. Turbidity currents have provided vital motivation to advance understanding of this class of flows because their enigmatic long run-out and driving mechanisms are not properly understood. Extant models assume that material transport by gravity currents is dynamically similar to fluvial flows. Here, empirical research from different types of particle-driven gravity currents is integrated with our experimental data, to show that material transport is fundamentally different from fluvial systems. Contrary to current theory, buoyancy production is shown to have a non-linear dependence on available flow power, indicating an underestimation of the total kinetic energy lost from the mean flow. A revised energy budget directly implies that the mixing efficiency of gravity currents is enhanced.

where h(x) denotes flow depth. It should be noted that it is assumed that flow parameters do not vary in time for quasi-23 equilibrium turbidity currents. Flow depth is defined as the height at which both flow velocity and concentration vanish, which 24 is 25 h(x) = max z| ⟨u⟩=0 , z| ⟨φ ⟩=0 , 26 where z| ⟨u⟩=0 is the flow depth above which flow velocity becomes negligible and z| ⟨φ ⟩=0 is the flow depth above which 27 concentration difference between the floa and the ambient fluid becomes negligible. In the ideal condition, z| ⟨u⟩=0 = z| ⟨φ ⟩=0 . 28 However, in the actual measurement of experiments, z| ⟨u⟩=0 is not always the same value as z| ⟨φ ⟩=0 for various reasons such as  (Table 1). 36 Thus, the accuracy of z| ⟨u⟩=0 and z| ⟨φ ⟩=0 varies between each experiment. In this study, the zero-velocity and zero-concentration 37 heights are carefully compared, and the more reliable one was employed as flow depth (see below for detailed procedures).

38
To describe the flow state, the following non-dimensional flow parameters are introduced. Firstly, the Richardson number, 39 Ri the ratio of the buoyancy and the flow shear stress, is given as, 40 Ri = −g(∂ ⟨ρ⟩/∂ z) where ρ is the flow density, ρ a is the density of the ambient fluid, and ∆ρ is the density difference between the flow and the 42 ambient fluid (⟨ρ⟩ − ρ a ). The shear velocity, u * , is calculated using the logarithmic law local to the bed, The skin drag coefficient, C D , is defined as 46 In this study, the particle settling velocity, w s , for relatively coarse material (> 100µm) is prescribed by an empirical formula where d 50 is the median grain size, and D s = (g∆ρ/ρν 2 ) 1/3 d 50 is the dimensionless particle diameter, ν is the fluid viscosity. 50 For the finer particles (< 100µm), Stokes' law is applied, where R is specific gravity and C 1 is a constant with a theoretical value of 18. Summary of compiled sources. The figures in the parentheses represent the number of measurement devices in the array. * 1 Numerical simulation of Reynolds-averaged Navier-Stokes model. * 2 Limited information about concentration profiles reported in this study. * 3 There are also experiments with sediment, but the availability of flow profiles is limited. * 4 Original reported values of particle size from hydrometer analysis. * 5 d 50 from Event 1, 4, and 5 from the original source. * 6 Velocity meters are special equipment that are designed for their particular study.
The compiled published sources are listed in Table 1. In total, 22 laboratory experiments, 1 numerical simulation, and 55 1 direct-observation source were gathered. Since each source uses different tools for vertical-profile measurement, the data 56 are compiled carefully. For example, the velocity measurement of flume experiments can be categorized into 5 different 57 types: micro-propeller current-meters, velocity meters, UVPs (Ultrasonic Velocity Profiler), ADVs (Acoustic Doppler Velocity 58 probes). For both UVPs and ADVs, each source used a different number of probes, sampling rate, and angle of mounting. 59 The vertical flow profiles reported in each source are based on interpolation and extrapolation, the methods differing between 60 the sources. To perform consistent analysis, here, only the original measurement points are extracted from each sources to 61 minimize the errors. Then, a consistent interpolation and extrapolation method are applied to recover the full vertical profiles.  Particle-size distribution analyses The potential error of particle size due to the different particle size analyses are carefully 65 considered. Settling velocity, w s , plays an important role in the required work done to keep sediment in suspension, B = gRΦhw s .

66
The particle size analysis methods in the compiled sources are either sieving combined with hydrometer method (SHM) or  Extensive comparison of SHM and LDM 28, 30-36 has shown a large discrepancy between the measured proportion of clay-size particles. Although the discrepancy becomes negligible for sand-size particles in most cases, when comparing the 74 size of the largest particles in sample of silt to clay-size particles, the error can be up to two orders of magnitude 33 . The 75 3/21 particle-size error of two orders of magnitude can result in up to 4 order of magnitude errors of settling velocity (Eq.(9)). Those 76 measurement discrepancies are attributed to the density difference and shape variation of particles 37 . 77 Supplementary Figure 1. a) The ratio of particle sizes as measured by the Laser Diffraction method to that of the hydrometer analysis method. Gray dotted line represents the fitted curve. b) Modified particle-size distribution of Tesaker 4 . In this study, particle size distribution from LDS is regarded as more reliable data than the SHM. Tesaker's experiment 4 78 used SHM to measure the particle size distribution of their clay material; the reported median particle size varies between 79 0.8-2.3 which is likely to be a significant underestimate. To mitigate this systematic error due to the difference of methodology, 80 in this study, the following modification of reported particle size is conducted.

81
Firstly, cumulative particle distribution curves of similar materials to that used by Tesaker (Kaolinite clay) are extracted 82 from the particle-size comparison studies between SHM and LDM 34, 35 . Since the particle-size discrepancy between SHM and 83 LDM could show different trends based on the type of clay minerals, only the datasets in which material is mainly composed of 84 Kaolinite were chosen. The cumulative particle distribution curve, where d i is the characteristic particle size of the i th bin, N is the totall number of bins, and f (d i ) denotes the fraction of the i th 87 particle size class which is given as where q(d i ) denotes the measured weight or volume of the particles in the i th size bin. Let r(d) denote the ratio of the cumulative 90 fraction of particles measured by LDM to the one measured by SHM, such that Supplementary Figure 1a shows the calculated best fit function r(d) from the compiled sources 34, 35 ; r(d) increases monotoni-93 cally as d increases until r reaches to unity at d ≃ 3 × 10 −2 .

94
Here, the relationship between r and log d for fine particles (d < 3 × 10 −2 ) is approximated by the linear curve fitting with 95 the least-square method (Gray dotted line in Supplementary Fig. 1a). Then, assuming r ≤ 1, the empirical formula of r is given   The first step is data formatting. Let P ⟨u⟩ = x ⟨u⟩,1 , x ⟨u⟩,2 , ...x ⟨u⟩,n and P ⟨φ ⟩ = x ⟨φ ⟩,1 , x ⟨φ ⟩,2 , ... respectively. x ⟨u⟩,i = (⟨u⟩| z=z u,i , z u,i ) and x ⟨φ ⟩,i = (⟨φ ⟩| z=z φ ,i , z φ ,i ) denote the i th coordinate of extracted points. The data points 117 are ordered from lowest to highest so that z i < z i+1 . Thus, z u,1 and z φ ,1 are the lowest original measurement heights and z u,n and 118 z φ ,m are the highest measurement heights. As boundary conditions, the following relationship is assumed for velocity profiles 119 in this study.
The next step is to estimate the flow height. The estimation of the flow height is performed in two steps: i) calculation 123 of z| ⟨u⟩=0 and z| ⟨φ ⟩=0 , and then ii) evaluation of flow height based on the comparison between z| ⟨u⟩=0 and z| ⟨φ ⟩=0 . Firstly, 124 interpolation of P ⟨u⟩ and P ⟨φ ⟩ is conducted using the PchipInterpolator function from Scipy, a Python package. This function 125 provides a piecewise cubic Hermite interpolation 38 . One of the advantages of this interpolation method is that it does not 126 overshoot even when the data is not smooth. Further, in this interpolation method, the original data points are always on the 127 interpolated curve, which means that potential oversimplification is minimal. Let f ⟨u⟩ (z) and f ⟨φ ⟩ (z) denote the interpolated 128 functions. Then, to obtain enough data points for processing, 500 vertically uniform data points are extracted from the 129 interpolated function, given as, where (a 1 , a 2 , a 3 ) and (b 1 , b 2 , b 3 ) are the coefficients of best-fit curves estimated by the least-square method. The conditions, 139 a 1 > 0 and b 1 > 0 are imposed to avoid unrealistic profile shapes. Then, z| ⟨u⟩=0 and z| ⟨φ ⟩=0 are obtained as, where φ ambient is the concentration of the ambient fluid. Then, the flow height is given by, 142 h = z| ⟨u⟩=0 when z| ⟨u⟩=0 ≤ z| ⟨φ ⟩=0 or z| ⟨φ ⟩=0 < z u,n , z| ⟨φ ⟩=0 when z u,n < z| ⟨φ ⟩=0 ≤ z| ⟨u⟩=0 .
143 z| ⟨φ ⟩=0 < z u,n sometimes happens when the number of data points in the density profile are inadequate. In this study, the original 144 velocity measurements are regarded as more reliable than the extrapolated concentration curve. Thus, when z| ⟨φ ⟩=0 < z u,n 145 ( Supplementary Fig. 2a,b), the extrapolation of the concentration profile is discarded and h = z| ⟨u⟩=0 is the flow depth. Then, 156 Finally, by conducting piecewise cubic Hermite interpolation against P ′ ⟨u⟩ and P ′ ⟨φ ⟩ , the full profile ⟨u⟩(z) and ⟨φ ⟩(z) were 157 obtained. The full profiles for each source obtained by the methodology described here are plotted in Supplementary Figure 3,  (TYPE III, see Table 1), the following empirical relationship 181 between Φ and the initial sediment concentration in the mixing tank, Φ 0 , are deduced, using orthogonal distance regression, 182 from TYPE II data where both parameters are available ( Supplementary Fig. 9). Then, Φ of TYPE III data is assumed to follow  Table 1 in the main text) are chosen to fill the gaps between the data points from the compiled sources ( Supplementary Fig.   197 11). A schematic and description of the experimental settings is given in the main text (Fig. 6 in the main text). Here, the 198 preparation of the sediment-water mixture, and the bed profile measurement which is conducted to measure the aggradation 199 rate, are elaborated.

200
Preparation of sediment-water mixture A large corn-shape mixing tank is built and used to generate turbidity currents (in 201 total 1.36 m 3 including the region below the mixer which cannot be used; the total effective volume for the experimental run   between total production, P loss and the log-law production,  Table. 1). The same methodology of data compilation for turbidity currents is applied for PDCs. For the entrainment rate  Supplementary Fig. 12a). The Yellow river data 46 are gathered from the figure from 57 , using a graph reading 239 software. Since the shear velocity values and raw velocity profiles of Yellow River data points are not available in the original 240 source, the average skin drag coefficient, C D of Yellow River data is estimated from the recent velocity measurement data 58 .

241
The reported flow parameters from each source such as flow velocity, volumetric flow concentration, and median particle size 242 are directly used for the calculation. For the detailed data compilation methodology for the rest of the data, see Dorrell et al. 59  of both explanatory and response variables, minimising the orthogonal distance between each data point and the fitted curve.

248
All explanatory variables in the linear regressions in this study are expected to exhibit error due to (for example) the limitation 249 of measurement tools such as UVP or siphon arrays. Thus, ODR is more suitable than OLS. In Table 1 in the main text, the 250 standard error of each estimated parameter is given, which is calculated from the square root of the corresponding diagonal 251 term in the parameter covariance matrix. To infer the goodness of fit of ODR, the coefficient of determination, R 2 is introduced.

252
R 2 is calculated from the sum of squares of residuals, RSS and the total sum of squares, TSS as,

254
where RSS and TSS for ODR is calculated as,

257
where (x i , y i ) denotes the coordinate of the i th observed data, d(x i , y i , f ) 2 denotes the squared orthogonal distance between 258 (x i , y i ) and the fitted curve f , and y denotes the mean value of y i .

259
Flow power plot To fit a curve to the data in figure 3a in the main text, we use the least squares fit for a parametric curve 260 fitting. To simplify notation, here we will use the notation x = log 10 (P f /Rghw s ), y = log 10 (Φ), so that the goal is to fit a curve 261 (x(t), y(t)) to the given set of data points (x i , y i ). The curves are written as polynomials in the parameter t, and we choose to 262 take x as a quadratic in t and y as a cubic in t, so that 263 x = 2 ∑ n=0 a n t n , y = This parametrization in the variable t has a symmetry of the form t → c 0 + c 1 t, which has corresponding transformations for the 265 coefficients a n , b n , that keeps the curve in the x-y plane unchanged. To remove this symmetry, we select two of the coefficients 266 by choosing t = 0 to correspond to the extremal value of x, a 1 = 0, and the rate of change of y around this point to be unity, 267 b 1 = 1.

268
To fit the curve to the data (x i , y i ), we recognise that we must not only determine the coefficients a n , b n , but also a set of 269 values t i which are the parameter values for the closest point on the curve to the data point. We do this employing the residuals
(31) 276 and we employ the method of least squares, seeking a local minimum of S(β β β ). Following the standard deviation for the method 277 of least squares, we can find a local minimum by iterating from one value, β β β k , to the next, β β β k+1 , using 278 (F F F T F F F + G G G T G G G)∆β β β = F F F T ε ε ε + G G G T δ δ δ (32) 279

284
2. Initialise a n and b m to random values between 0 and 1, except for a 1 = 0 and b 1 = 1.